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Abstract 


An efficient numerical approach for the design of optimal aerodynamic shapes is presented in this paper. 
The objective of any optimization problem is to find the optimum of a cost function subject to a certain 
state equation (governing equation of the flow field) and certain side constraints. As in classical optimal 
control methods, the present approach introduces a costate variable (Lagrange multiplier) to evaluate the 
gradient of the cost function. High efficiency in reaching the optimum solution is achieved by using a 
multigrid technique and updating the shape in a hierarchical manner such that smooth (low-frequency) 
changes are done separately from high-frequency changes. Thus, the design variables are changed on a 
grid where their changes produce nonsmooth (high-frequency) perturbations that can be damped efficiently 
by the multigrid. The cost of solving the optimization problem is approximately two to three times the 
cost of the equivalent analysis problem. 


Nomenclature 

F cost function 

fk fcth shape function 

.1/ free-stream Mach number 

n unit normal 

t unit tangent 

TJ y free-stream velocity 

y v : L ^-coordinate of the upper and lower surface of the 
airfoil 

<v amplitude of shape functions (design variables) 
o direction of change of a 

q ,V l components of o (upper and lower surface 
amplitudes of the k th shape function) 

T circulation 

7 ratio of specific heats 

e magnitude of change of o 

C angle of attack 

0 6 corrected for Mach number 
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6 angular position of a far-field location 

A Lagrange multiplier 

fik fcth component of the gradient of F 

p density 

o full velocity potential 

o ,q target potential 

T coefficient of the delta function 

I. Introduction 


The aerodynamic designer seeks an airplane shape that 
achieves the best aerodynamic performance while taking into 
account the trade-offs with other disciplines, namely, struc- 
tures, propulsion, stability, and control. Although significant 
progress has been made in developing computational meth- 
ods for fluid flow analysis, the methods for the design and 
optimization of aircraft configuration are still in their infancy. 
Among the many techniques that have been developed, the 
inverse design method 1-1 ’ is perhaps the most widely used. 
This method is severely restrictive because its depends on 
the experience and knowledge of the designer to establish 
desirable velocity or pressure distributions. In addition, the 
method does not lend itself to the imposition of constraints. 
The next most popular approach is what we refer to as the 
“loosely coupled optimization (LCO).” In this approach, the 
analysis problem is solved many times to find the gradient 
of the cost function with finite differences. This gradient is 
used by a “black box” optimizer to find the best change to the 
design variables. The computational cost of this brute-force 
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method is usually exorbitant. Examples of this approach are 
found in Refs. 7-10. The LCO method can be improved 
by analytically evaluating the sensitivity derivatives that are 
needed to update the design variables. 11 Usually, this re- 
quires the inversion of an extremely large matrix. For three- 
dimensional problems, the size of this matrix can render the 
method impractical with current computer technology. In the 
past several years, the adjoint variable method has gained 
popularity in the area of aerodynamic design. This method, 
which we refer to as “tightly coupled optimization (TCO),” 
requires the solution of an adjoint problem, in addition to the 
analysis problem, to evaluate the gradient of the cost func- 
tion. The complexity of the adjoint problem is equivalent to 
that of the analysis problem. However, the adjoint problem, 
along with the analysis problem, must be solved many times 
to reach convergence. This approach has been discussed in 
Refs. 12 and 13. Although the efficiency of TCO is much 
higher than that of LCO, this procedure can also become 
prohibitively expensive for practical aerodynamic design and 
optimization problems. 

The One-Shot method 1415 overcomes the unacceptable 
cost of the existing design and optimization procedures. It 
brings the cost of design and optimization to the same order as 
that of a single analysis. High performance is achieved both 
by exploiting the property of the partial differential equations 
(associated with the scales (frequency) of the errors) which 
govern the physics of the flow and by the efficient damping 
of high-frequency error components with multigrid. Consider 
the subsonic flow over an airfoil profile. The change in the 
shape of the profile of a given wavelength produces changes 
of the same wavelength in the solution. These changes 
penetrate into the flow field only up to a distance that is 
proportional to the wavelength of the perturbation. Thus, 
while the high-frequency changes in the shape of the airfoil 
produce changes in the solution that are of high frequency 
and remain local to the neighborhood of the airfoil, the 
smooth (low-frequency) changes in the shape produce smooth 
changes to the solution and are global in nature. Typically, 
any relaxation scheme quickly damps the high-frequency 
components of the error on a grid. Multigrid efficiently 
damps the whole spectrum of error components by relaxing 
the governing equations on a sequence of grids of varying 
resolution. 

Therefore, the basic idea of the One-Shot method is to 
change the shape of the airfoil profile in a hierarchical manner 
such that smooth changes are made separately from high- 
frequency changes. Because each of these changes involves a 
different scale, the governing equation of the flow field can be 
solved efficiently on grids of appropriate resolution. Thus, the 
flow field due to smooth changes in the shape of the airfoil is 
solved on coarse grids, and the flow field due to increasingly 
high-frequency shape changes is solved on increasingly fine 
grids. This breaks the optimization procedure into a sequence 
of suboptimization problems, each of a given scale; therefore, 
the problem is well conditioned. The resulting optimization 
procedure is very efficient because the work on a particular 
scale is done on the appropriate grid. (Ill conditioning results 


from working on many scales simultaneously.) The One- 
Shot method is implemented within a full approximation 
scheme (FAS) full multigrid (FMG) algorithm. The solution 
process starts on the coarsest grid, where only the smooth 
component of the shape function is updated. This solution 
is interpolated to the next finest grid, where it serves as 
an initial approximation of the solution on that grid. This 
process is continued until the finest grid is reached. Thus, 
smooth (low-frequency) shapes are updated on coarse grids; 
high-frequency shapes are updated on finer grids. The fine- to 
coarse-grid transfers are designed such that the optimization 
problem at each grid level is driven by the fine-grid residual. 
The resulting algorithm has an estimated overall cost that 
ranges from two to three times the cost of the analysis 
problem. 

The successful application of the One-Shot method to 
the aerodynamic shape design problem was first reported in 
Ref. 15. The capability of the method was demonstrated by 
using the small-disturbance potential equation as the govern- 
ing equation of the flow field. However, in that study, the 
issue of updating the grid was avoided. In the present study, 
the full potential equation is used as the governing equation; 
hence, the grid must be updated as the shape changes. In this 
work, the adjoint equation and the corresponding gradient of 
the cost function are derived. The solution procedure and 
some typical results are also presented. 

II. Design of Optimal Airfoil Shapes 


The design of optimal airfoil shapes is a constrained 
minimization problem. The objective is to find the optimal 
shape of the airfoil that will minimize a cost function F 
subject to the state equation of the flow field and the side 
constraints. 


The State Equations 


The analysis problem, defined by the state equation, 
consists of finding the flow over a specified shape for a given 
free-stream Mach number and angle of attack. In order to fo- 
cus on the optimization procedure, the flow model considered 
is the subsonic potential flow over an airfoil profile. 

Consider the steady irrotational flow past a two- 
dimensional airfoil. 1617 The governing equation of the flow 
field, known as the full potential equation, is 

div(pV</>) = 0 (1) 

The boundary condition on the airfoil is 

V< f> n = 0 (2) 

At infinity, the boundary condition is 

V<t> = Uoc (3) 
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Fig. 1. Computational domain. 

For the Kutta condition, the circulation I around the airfoil 
is such that 


where t.e. refers to the trailing edge of the airfoil. To satisfy 
mass conservation across the cut, derivatives of the potential 
normal to the cut are required to be continuous. 

At the far-field boundary, the circulation modifies the 
velocity as follows: 


r 

V<j 4 • n = Uc c, • n H V0 • n (9) 

27T 

where 

0 = 27r — tan -1 (a/1 — M tan d'j (10) 

and 8 is the angular position of a far-field point. For conve- 
nience, n is the unit normal on the boundary. The far-field 
boundary condition given by (9) is consistent with the infin- 
ity condition stated by (3). 


The Design Variables 


The airfoil is represented as follows: 


the velocity at the trailing edge is finite and continuous 

(4) 

In these equations. <f> = <f>(x,y) is the full velocity potential, 
p = p((f>) is the density, n is the unit normal, and U 1>Cl is the 
free-stream velocity. The density p is given by 

1 

7-1 

(5) 


P = 


1 - 


JV*| 2 - 1 ) 


where Mo% is the free-stream Mach number and 7 is the ratio 
of specific heats. If ( is the angle of attack of the airfoil, then 
the free-stream velocity is given by 


Uc c = Uco [cos( C)*4- sin( ( )j] (6) 


where i and j are the unit vectors in the x and y directions, 
respectively. 


The Computational Domain 

The computational domain is shown in Fig. 1. The 
interior of the flow field is denoted by Q ; the upper and lower 
surfaces of the airfoil are denoted by U and L, respectively. 
The far-field boundary, located at a finite distance from the 
airfoil (30 to 50 airfoil chord lengths) is denoted by O. To 
impose the Kutta condition around the airfoil, an artificial 
boundary or cut that begins at the airfoil and extends to the 
far field is introduced. A jump in potential that is equal to T 
is allowed across the cut. For convenience, this cut is chosen 
to emanate from the trailing edge of the airfoil. The top and 
bottom sides of the cut are denoted by T and B. respectively. 
The jump across the cut can be written as 

<t> T -<t> B = T (7) 

The value of the T is determined by requiring that the 
velocity perpendicular to the trailing edge bisector be equal to 
0 at the trailing edge. A good approximation for T is given by 

r = - 4>t.e. (8) 


K 

y V = y ^akfk(x) 

k = 1 

(0 < X < 1) (11) 

K 

y L = alfk(x) 

k— 1 

where o): T and a\ are the amplitudes of the shape functions fk 
on the upper and lower surfaces of the airfoil, respectively. 
The design variables ctk must be determined to obtain the 
optimal shape of the airfoil. Let a denote a vector whose 
elements are the design variables. That is, 

r U U U L L l] T 

n := ,«2 , ( 12 ) 

The functionality of the shape functions will be presented 
later. 


The Optimization Problem 

The model problem chosen is the design of an airfoil 
shape that can match a given target potential. Given a target 
potential distribution o r, around an airfoil, the objective is to 
find a that will minimize 

F[a, <t>(a)] = j (<t> — <t>o) 2 da (13) 

U + L 

subject to the state equations, where da (which is an element 
of the airfoil) can be written as 

da 2 = dx 2 + dy 2 (14) 

Note that the choice of this particular cost function does not 
make it an inverse-design problem. Unlike inverse-design 
problems, the minimization is done over a finite number of 
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design variables. This approach also can be used, for ex- 
ample, to find the optimal shape of an airfoil that has the 
minimum D/ L (drag/lift) subject to geometric and aerody- 
namic constraints. 

To make the presentation of the derivation of the ad- 
joint equations simple and easy to understand, the flow is 
assumed to be incompressible (i.e., .1/ = 0); therefore, 

p = 1 . In this case, the full potential equation reduces to the 
Laplace equation. Also, no side constraints are considered in 
this derivation. Therefore, the specific optimization problem 
considered here is 


mm 

OL,<fi 


j (<t> ~ <t>o 


" da 


U + L 


(15) 


subject to 


The Step Size. The step size e is determined by a line 
search. The objective of the line search is to find e such that 
| V a F (at + e«, <j> + e<j>] | is a minimum. That is, 


c> || V n F ( at + at, <j> + e 

d£ 


= 0 


( 21 ) 


By expanding in Taylor’s series, taking the derivative with 
respect to e, and neglecting the 0(e 2 ) terms, the step size 
becomes 


[V n F( n ,<t>)] T VlF( n ,<t>)a 
a T [ V l F ( a, 4> )] T V l F ( a, 4> ) a 


( 22 ) 


di v(V<t>) = 

0 


in Q 

(16a) 

<\ 

-e- 

S 

II 

0 


on the airfoil 

(16b) 

<1 

-e- 

S 

II 

Ucc 

r 

■n- 1 V0-n 

2tt 

in the far field 

(16c) 

1 

-Sk- 

td 

II 

r 

along the cut 

(16d) 


where V^.F is a symmetric matrix and is often referred 
to as the Hessian. Note that V^F includes the variation 
with respect to Computation of the Hessian is expensive; 
the cost is proportional to the number of design variables. 
However, V;, / A can be evaluated relatively easily with 
finite differences as follows: 


where T is given by (8). 


The Minimization Process 

At some initial a, any minimization process seeks to 
find a descent direction at and a step size e in which to 
change at such that 

F (a + eat, <f> + etfrj < F(a, <f>) (17) 


V 2 L(a, 4>)ol = 


V a F (a. + £Ot , tf> + £<f> ) — V,-»L(rt, i j>) 


(23) 

where e is a trial perturbation. To find the step size, 
the design variables are first perturbed with an arbitrar- 
ily small e, and the new values of the state variables that 
satisfy the state equations are determined. Next, the new 
gradient V a F (at + £Ct :, <j> + e<j>\ is evaluated, followed by 


V 2 _F(a, <f>)a. Then, the step size is determined with (22). 


where £tj> is the corresponding change in 4 that satisfies the 
state equations. This process is repeated several times until 
a minimum is reached. 

The Descent Direction. A descent direction at can be 
determined as follows. The Taylor series expansion of F 
about at and <j> can be written as 


F (a + eat, <f> + e<j>^ = F(at, <t>) + £C\ T V a F(n, tf > ) + 0(e 2 ) 

(18) 


where 




Equation (18) clearly shows that if 

- _ V^(a, 4 >) 


\V a F(a,4>)\ 


(19) 


( 20 ) 


The Adjoint Equations 

As stated earlier, the objective of the optimization pro- 
cedure is to seek a descent direction and a step size in which 
the design variables can be changed so that the cost func- 
tion is decreased. To determine the descent direction and the 
step size, the gradient of the cost function with respect to the 
design variables V d f (given by (19)) must be evaluated. Ef- 
ficient and accurate evaluation of the gradient of F is one of 
the important but difficult steps in any minimization scheme. 
Equation (19) is not very useful for determining the gradi- 
ent of F because the efficient and accurate determination of 
d<f>/dot generally is difficult. However, the adjoint method 
offers an elegant means for evaluating the gradient easily 
and accurately. The derivation of the adjoint equations is 
presented below. 

Let the design variables be perturbed such that 


then (17) is satisfied. Equality occurs in (17) at the minimum 
when V a F{a* , <f>*) = 0, where o* is the optimum value of 
the design variables and <t> * is the corresponding value of the 
state variables that satisfies the state equations. Therefore, 
to obtain the descent direction, the gradient of F must be 
evaluated. 


O' — « + e« (24) 

where £at is the change in at\ £ and at are the step size 
(magnitude) and direction, respectively, of the change in o.. 
Figure 2 shows the domain after the perturbation, where U 
and L denote the upper and lower surfaces, respectively, of 
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Fig. 2. Domain after perturbation. 

the new airfoil and Cl denotes the new domain. The shape 
of the resulting airfoil t/ u,L and the corresponding potential 
a that satisfy the governing equation and its boundary con- 
ditions in the new domain can be written as 


If we introduce a Lagrange multiplier A and use (29a), 
then (27) can be written as 


a T V n F 


2{<f> — <f>o)<f> y yda + I (cj> — cj> of dcr 


U + L U + L 

+ 2( cf> — <t>o )<j>da + jj div ( V<t> ] XdFl 

u+L n 


l + m 


(31) 


If we integrate by parts, the last integral can be written as 


.// 


cliv ( VS ) XdQ, 


II 


div( VA)+/fi 


- I a(v^ • njdr + j (V\-n,)<j>dr 

T T 

(32) 

where the unit normal n, points into the flow field fi; dr is 
an element on r , which is the path of integration around the 
domain Q and can be expressed as 


_U,L 

y 


U,L 

= y 


+ sy 


(25) 


7- = L + U + T + 0 + B (33) 


and 


4> = 4> + e4> 


If the integrals along r are split into different components 
(26) and substituted into (31). then we can write 


where ey represents the change in the airfoil shape and s a 
represents the corresponding change in the potential. We can 
show from (13) that 


a T V n F 

= I 2(4 

U + L 


4>o)4>yyd(j 


+ j (<t> 

U + L 


<t>of 


y y 

1 + yl 


da 


+ j 2 (<t> — <t>o)<j>da 
u+L 

(27) 

where y , = dy/dx. The objective of this approach is to 
eliminate those terms that have <j>, where 


1 d<t> . 
<t> = n 

da 


(28) 


a T V n F 

= I 2( <t> 

U + L 


4>o)4> v yda 


+ j (<t> 

U + L 


4>o f 


y y 

1 + yl 


da 


F j " 2 (<t> — <f>o)<f>da +- 

ff div(VA )<f,dU 

U + L 

a 


- 1 x(^V<j>-njda 

+ 1 (VA 

• n)<f>da 

U + L 

U + L 


- f x(v<f>-nfjdr 

+ / (VA 

■ n)<t>dr 

T + B 

T + B 


- A (^Vcf> ■ nfjdr 

+ /.VC, 

n)<f>dr 


o o 


(34) 


From the governing equation and its boundary conditions 
(16), we can show that 


div^V^ = 0 

in Q 

(29a) 

V+n = (yV<t>-t) x 

on the airfoil 

(29b) 

f 

Vd • n = — VO-n 

')t 

at the far field 

(29c) 

1 

td 

II 

^ t 

along the cut 

(29d) 

where 

f = ^L. 

tb 

— 4>t.e. 

(30) 

and t is the unit tangent. 




Because V<t> is continuous across the cut and n points in 
opposite directions along the top and bottom boundaries of 
the cut, we can write 

j \(v<f>-n)dr = j (\ T - njdr (35) 

T + B Cut 

If we assume that VA is continuous across the cut, then we 
can write 

j " (V A • n)<f>dr = f VX-ndr (36) 

T + B Cut 
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If we use (35), (36), and (29b-d), then equation (34) can be Because div(VA) = 0 in Q, we obtain the following 

written as from the divergence theorem: 


a T V n F 


2 (4> — 4>o)4> v yda + 


U + L 


I [<f> — <t>o) 2 77 da 


U + L 


1 + ys 


j 2( <f> — <t>o )<f>da + jj div(VA)+(Q 

+ L £1 

J ■ t) x da + j (VA • n)<t>d<T 


U + L 


U + L 


j (a t - A b ) (v^ • njdr + f j V\- ndr 


- — / A( V0 • n)d.T + j (VA • n)<t>d,T 
o o 

(37) 

If we substitute for f from (30) and rearrange, then (37) 
becomes 


n T V n F 


j 2(<f> -<t> 0 )4> v yd<j + I (<t> 

+ L U + L 

j \{ yV<f> ■ t)J.a + jj div(VA )<f>dU 


i \ 2 VxVx , 

— 0o J — — -da 


1 +y£ 


U + L 


+ / [VA • n + 2(0 — 0o)]0^o 


U + L 


U TL 
t.e. ~ <Pt.e 


^ VA • ndr — j AV0 • ndr 


+ j (VA -n)4>dr- J (a t - A B )(v+ti)t/7- 

O Cut 

We choose A such that 


(38) 


div(VA) = 0 in Q 

VA • n + 2((j> — <f>o) — T6(i — it, e .) = 0 on L 

VA • n + 2((j> — <f>o) + T6(x — xt.e.) = 0 on U 

VA • n = 0 in the far field 

A T — A B = 0 along the cut 
(39) 


where 


j VA • ndr ~ j - 


T = / VA • ndr / AV0 • ndr (40) 

27T 


and 5 denotes the Dirac delta function. Equations (39) are 
the adjoint equation and its boundary conditions (also called 
the costate equations). These equations are similar to the 
linearized state equations. The size of the system is the same 
as the size of the state equations and can be solved with the 
same technique used to solve the state equations. 


j VA • ndr = 0 


( 41 ) 


Therefore, so that (39) has a solution, we can show that 

j j ■ ■ <t>o)d(j = 0 (42) 

U + L 

Equation (16) clearly shows that a constant can be added to 
n. We can choose this constant o , such that 


j (<t> + 4 > c — <t>o)d(j = 0 

U + L 


Therefore, 

[(</> — <t>o)d<r 

A U + L 

= fj , ; — 

U + L 


(43) 


(44) 


The Gradient of F 

If (39) is substituted into (38), then it reduces to 

& T V n F 


2(<t> — <t>o)<t> y yd<r + / (<t> — <t> o 


U + L 


U + L 


\ 2 VxVx 

i+yl 


-da 


— j \{yV<f> ■ t) x da 

U + L 

If we integrate the last integral by parts, then we get 


(45) 


A (yV4> ■ t) x da 


U + L 


j yV<t> ■ t\ x da + j 


yV tj> ■ t\ x da + I yV <j> ■ tX ( l a 

1 + Vx 


(46) 


U+L U+L 

If (46) is substituted into (45), then we can write 


a T V a F 


$_(<!> — <t>o)<t> v yda + / (<t> — 4>o ) 2 nxflx n da 


U + L 


U + L 


l + ys 


j yV <f> ■ t\xda + j 


+ I yV tj> ■ tXxda + / yVcj) ■ /A 11 11 da 


U + L 


U + L 


1 + yi 


(47) 

If we substitute for y from (11), then (47) can be written as 

I\ K 

n T V a F = +Y, a ^d‘i ( 48 ) 
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where 


/4 T,L = f (<t> -4>o) 2 1 + y 2 (fk) x d<r 

U,L 


+ 



' t\x + S7<t> ' tX 


y xVxx 

1 + 2/1 


fkda 


(49) 

Equations (49) are the components of the gradient of F. 
When o satisfies the state equations (16) and A satisfies the 
costate equations (39), then the components of the gradient 
of F can be evaluated with (49). Because V„F = 0 at the 
minimum, we can clearly see that 


hk = 0 1 

Uk = 0 J 


for k = 1,2, ...K 


(50) 


A Design Strategy 

Figure 3 shows a typical design strategy. In this process, 
at some initial conditions the state and adjoint equations are 
solved, and the gradient of F is computed. If the gradient is 
equal to 0, then a minimum has been reached and the iteration 
is terminated; otherwise, the new descent direction a and 



Fig. 3. A design strategy flowchart. 


the step size e are computed, and the design variables are 
updated. The iteration is repeated until the gradient vanishes. 
The cost of this strategy can be estimated as follows. Let 
the cost of solving the state equations be equal to K. The 
cost of solving the adjoint equation is at most equal to K. 
Let the number of design iterations required be N. Therefore, 
the total cost of doing the optimal design is approximately 
2KN with N, at best, equal to the number of design variables. 
In practice, especially for nonlinear problems, N is many 
times the number of design variables. A factor of 100 is not 
unrealistic. One way to bring the total design cost down is 
to reduce the magnitude of K. One of the most practical and 
proven ways of achieving this is by using multigrid. Here, 
a multigrid scheme is used to relax the state and adjoint 
equations. At the end of one or several multigrid cycles, 
the optimizer is called and the design variables are updated. 
In this process, the design variables are updated only on the 
finest grid. A schematic of this strategy is shown in Fig. 4. 



Fig. 4. A multigrid strategy. 
III. The One-Shot Method 


The One-Shot method goes one step further by embed- 
ding the design process within the multigrid cycles. This 
method essentially makes N = 1. Thus, the cost of optimal 
design is approximately equal to 2 K. In this method, high 
efficiency is obtained by exploiting two key phenomena: the 
ability of multigrid to efficiently reduce high-frequency com- 
ponents of the error due to a perturbation and the nature of 
propagation of perturbations in a flow field. These phenom- 
ena are explained below. 


Multigrid Efficiency 


In any relaxation (smoothing) process, the high- 
frequency error components of the space discretization op- 
erator of the differential equation under consideration are 
generally damped in a few iterations. The low-frequency 
components are the slowest to be damped. Consider a one- 
dimensional domain of length L discretized into N cells of 
uniform grid spacing h = L/N , where the grid index ranges 
from 0 to N . This grid will be referred to as the h grid. If 
we assume periodic boundary conditions, then the error at 
the nth grid point can be written in Fourier series as 

N 

e„ = A 3 e ' e,n (5!) 

j = -N 
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where Aj is the amplitude of the jth harmonic and i = \/— T. 
The phase angle 6 can be written as 


change in shape to /+ /, the governing equation for change 
in potential <f> is 


The phase angle covers the domain (— 7r,7r) in increments 
of 7 t/N. The value |#| = 7r corresponds to the highest 
frequency that is visible on this grid, namely the frequency 
of wavelength 2 h . If a coarse grid ( H grid) is constructed by 
removing every other grid point of the h grid, then the highest 
frequency that is visible on this grid corresponds to |6* | = 7r/2 
(i.e., the frequency of wavelength 4 h = 2 H). Therefore, the 
frequencies that correspond t0 7r/2<|#|<7r and are visible 
on the h grid cannot be represented on the H grid. These 
frequencies are considered to be high frequencies on the h 
grid, and the relaxation scheme can damp these frequencies in 
a few iterations. The remaining frequencies in the spectrum, 
which correspond toO < |#| < 7r/2 and are well represented 
on the H grid, are referred to as low frequencies on the li 
grid. The frequencies that are visible on the H grid can also 
be separated into high and low frequencies, based on how 
well they are represented by the next coarsest grid. The high 
frequencies that correspond to the H grid can be damped 
quickly by a few iterations of the relaxation scheme on this 
grid. 

In the multigrid method, 1718 high efficiency is obtained 
by relaxing the discretized equation on successively coarser 
grids, where the high-frequency error components that cor- 
respond to each grid are damped efficiently. In the design 
process, high efficiency is obtained by changing only those 
design variables that produce high-frequency perturbations in 
the flow field on any grid. Therefore, the basic premise of the 
One Shot method, on any grid, is to make changes in the de- 
sign variables that produce high-frequency perturbations in 
the flow field. 


The Effect of Airfoil Perturbation on the Flow Field 


The other phenomenon that is exploited by the One- 
Shot method has to do with the way in which a disturbance 
is propagated in a flow field. In a subsonic flow, for example, 
a smooth perturbation is propagated through the entire flow 
field and a high-frequency perturbation is felt only in a small 
neighborhood around the source of the perturbation. That is, 
high-frequency components of the perturbation decay rapidly 
away from the source. This phenomenon is illustrated in the 
following analysis. 

Consider the small-disturbance potential equation in the 
half-space 0 < y < oo, — oo < # < oo. If the flow is 
incompressible, the governing equation is 

V 2 <f> = 0 (53) 


and the boundary condition applied at y = 0 is 

d<t> _ df 

dy dx 


(54) 


where f(x) is the shape of the boundary over which the 
flow must be determined. If a + o is the potential due to a 


v 2 4 > = o 


(55) 


and the boundary condition at // = 0 is 


d± = d£ 

dy dx 


(56) 


Let 

df _ 
dx 


(57) 


where ui is the frequency of the perturbation. A solution 
to the governing equation (55) that satisfies the boundary 
condition is 


= e“ H V ; 


(58) 


The magnitude of <f> is 

|^| = e“ Hs ' (59) 

Figure 5, which is the plot of (59) for a few select frequencies. 



Fig. 5. o versus y. 


shows that the region where <fx is large becomes thinner as 
the frequency increases. Let y* be a location where o is less 
than some small it. That is. 


4>(io,y*) <£ 


„-/3 


(60) 


If we substitute for o, then 


Therefore, 


e -Ms/* =e -/3 

(61) 

» — In ( e ) 

M 

(62) 


Equation (62) clearly shows that as the frequency of the 
perturbation u> increases y* decreases. Table 1 shows y* 




for a few select frequencies when e = 10 4 . For the discrete 


The Shape Functions 


Table 1. y* Versus w 



problem, (59) can be written as 

|^| = e -(h|hK»/h) =e -|e|U-i) {j = f,2,...J + 1) 

(63) 

where itf J < 8 < it is the frequency scaled to the grid spac- 
ing h. Figure 6 shows the response to different frequencies 
for the discrete problem. Table 2 shows the grid location j * , 
beyond which |^| < 10 _4 . It shows that the high-frequency 
perturbations are significantly damped by about the fifth grid 
point (j = 0 is the first grid point). 
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Fig. 6. o versus j. 


As presented earlier (section 3), the airfoil is represented 
as follows: 

y V ’ h = '^2 a k’ h fk(x) (64) 

k— 1 


where o): T and a\ are the design variables and fu are the 
shape functions. As explained in the previous two sections, to 
obtain high design efficiency, the changes in the design vari- 
ables on a grid should produce nonsmooth (high-frequency) 
perturbations in the flow field. This is achieved by using a 
set of orthonormal functions as shape functions. Orthonor- 
mal functions are increasingly oscillatory. Each of them is 
assigned to a grid where a change in the amplitudes causes 
nonsmooth perturbations in the flow field. Often, basis func- 
tions that correspond to some known airfoil shape must be 
used. If these functions are not orthonormal, then the cor- 
responding orthonormal functions can be determined by a 
Gram-Schmidt process. A Gram-Schmidt procedure for or- 
thonormalization can be developed with the property of or- 
thonormal functions, namely. 


j fm ( X I./ /, f X ) dx — 0 (in 
o 

l 

j f m ( X ) — 1 


Let yk(x) be the functions that are not orthonormal. First, the 
orthogonal set //, ( j ) is found from the following relations: 


fi(x) = gi(x) 

fcix) = (J 2 (x) + a. 21 .fl(x) 


k-1 

fk(x) = g k (x) + E (t k m f m f 1 ) 

777 = 1 


( 66 ) 


Table 2. j* Versus 9 



In the One-Shot method, a shape function is perturbed 
on a grid where it produces high-frequency error compo- 
nents. These errors penetrate only a small distance into the 
flow field. Hence, they can be quickly damped by a few re- 
laxations of the discrete equations in a small neighborhood 
around the airfoil. 


where 

i 

/ gk(x)f m (x)dx 

a km = — - — j (67) 

f fm(x )dx 
0 

Finally, the orthonormal functions are found by normalizing 
fk(x) as follows: 



The Gram-Schmidt process described above can be pro- 
grammed in symbolic language to find the expressions for 





fk, or it can be implemented by numerical integration, in 
which case the shape functions are defined as an array of 
numbers. 

As an example, consider the NACA 0012 airfoil, which 
is defined by 


y v = T: Ihgkjx) (0 < x < 1) 

k=i 



(69) 


where A and gu are given in Table 3. The NACA 0012 
shape has been slightly modified to ensure that it closes at 
the trailing edge. The same shape can be expressed in terms 
of the orthonormal functions as 

4 

y V = T akfk(x) (0 < x < 1 ) 

hi " ' ( 70 ) 


where the orthonormal functions fk of the basis functions 
and their corresponding amplitudes n /, are given in Table 4. 
The orthonormal shape functions are shown in Fig. 7. Note 
that the number of zeros in fk is equal to k + 1 . 


Table 3. Shape Functions and 
Amplitudes of NACA 0012 


k 

0 k 

Qk 

1 

0.17814 

\fx — X 

2 

0.10128 

i(l — X ) 

3 

-0.10968 

IN) 

i 

4 

0.06090 

S3 

CO 

1 

S3 


Table 4. Orthonormal Shape Functions 
and Amplitudes of NACA 0012 


k 

a k x 10 4 

fk 

1 

439.474 

5.47723.9! 

2 

28.2339 

14.7573(92 - 0.928571.9i) 

3 

-5.85699 

54.7884(93 - 0.901236.92 
+ 0.432099.91) 

4 

2.85283 

213.472(94 - 1.27406.93 
+0.50401192 - 0.164439.91 ) 



Fig. 7. Orthonormal shape 
functions of NACA 0012 airfoil. 


The One-Shot Design Strategy 


In the One-Shot method, the optimizer is embedded 
within the multigrid cycle as shown in Fig. 8. The de- 
sign variables are updated on a level where the correspond- 
ing shape functions produce high-frequency error compo- 
nents. In general, the low-frequency shape functions are 
updated on coarse levels, and higher frequency functions 
are updated on finer grids. For example, the design vari- 
ables o? and a\ are updated on the coarsest grid 8 h\ 

U U U U L L L , L , . , 

a i , a 2 , a 3 , a 4 , a 1 , a 2 , a 3 , and 0 4 are updated on 
the next finest grid 4 h. Some overlap of the design vari- 
ables is permitted. Thus, a? and a \ are updated on grid 4 h 
also. None of the design variables are updated on the finest 
grid h. The cost of solving the state or the adjoint equations 
on a coarse grid is only one-fourth of the cost of solving 
them on the next finest grid. Because the shape functions are 
perturbed only on levels where they generate high-frequency 
errors, a local relaxation around the airfoil is sufficient to 
damp out the errors, which reduces computing costs. There- 
fore, the overall cost of the design is dominated by the cost 
required to solve the state and adjoint equations on the finest 
grid. The total cost of the design process is approximately 
two to three times that of one analysis. 



Fig. 8. The One-Shot strategy. 
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The Discretization and Solution Procedure 

The State Equations. The computational domain is 
discretized with an O type of grid. The governing equation 
and its boundary conditions cast in curvilinear coordinates 
are discretized with the finite-volume approach. The Gauss- 
Seidel line-relaxation scheme is used to form the tridiagonal 
systems of equations in both curvilinear coordinate directions. 
These systems are solved with the Thomas algorithm. Note 
that the tridiagonal system is periodic in the direction that 
is around the airfoil. A FAS multigrid scheme is used to 
accelerate the convergence rate of the solution. The FMG 
process is used to obtain a good initial solution on the finest 
grid. 

The Adjoint Equations. The adjoint equations are dis- 
cretized and solved in the same manner as the state equations. 
As in the case of the state equations, a FAS multigrid scheme 
and the FMG process are used to accelerate the convergence 
rate of the solution. 

The Gradient of F. The gradient of the cost function 
involves only quantities on the airfoil. These quantities are 
discretized in a manner that is consistent with the discretiza- 
tion of the state and adjoint equations. The gradient is trans- 
ferred to the coarse grid in a FAS manner. 

Updating the Grid. During the design process, the 
grid is updated by moving only the grid points close to 
the airfoil and linearly decaying the change at the airfoil 
in this neighborhood. The outer boundary of this region is 
determined as follows. Let 

|/max — ijmax(ej) (73) 

where t] is an arbitrary constant; i] = 10 in this study. Among 
the grid lines that go around the airfoil, the one that is nearest 
to the y max location is taken to be the outer boundary of the 
region within which the grids are changed. The entire grid is 
regenerated at the beginning of each FMG stage also. With 
this approach, by the time the FMG process reaches the finest 
grid, only a few lines around the airfoil must be moved. 

IV. THE RESULTS 

Test Case 1. 

As our first test problem, we recover the NACA 0012 
airfoil shape using the potential distribution obtained from 
the analysis of NACA 0012 at an angle of attack of 0° 
and .1/ = 0 as the target potential <f>o- Figure 9 shows 

the computed C p distribution obtained from the analysis run. 
A five-level W-cycle multigrid with 128 x 64 cells on the 
finest grid was used. The FMG process was used to obtain 
a good initial approximation for the finest grid. The analysis 
converged to machine zero (< 1 0 — 1 0 ) in 10 multigrid cycles. 

The design run was similar to the analysis run. During 
the design process, both the state and costate equations were 
relaxed at any multigrid level. The shape functions used were 
the orthonormal functions based on the NACA 0012 shape 
functions. The design variables were distributed such that on 




Fig. 9. Computed C p distribution for NACA 0012. 

the coarsest level (8 x 4) only of and of were updated. 
On the next finest level (16 x 8), all the design variables 
(<> ) were updated. None of the design variables were 

updated on the next three levels, including the finest level. 
Thus, most of the design overhead was limited to the two 
coarsest grids. The FMG process was used to obtain a good 
initial approximation of the shape for the finest grid. Figure 
10 shows the results of this run. The residuals of the state and 
costate equations and the gradient of the cost function reached 
machine zero in 12 multigrid cycles. The cost function at 
convergence was equal to 3 x 10~ 13 , which indicates that 
NACA 0012 was indeed recovered. 



0.0 0.2 0.4 0.6 0.8 1.0 


X 

Fig. 10. Test case 1. 

Test Case 2. 

For test case 2, we selected the airfoil FX 60-126/1, 
a cambered airfoil whose coordinates are tabulated in Ref. 
19. Figure 11 shows the C p distribution for this airfoil at 
an angle of attack of 0° and M = 0. This airfoil is not 
smooth, which is reflected in the computed C p distribution. 
Using this solution as the target, we tried to recover the shape 
with the NACA 0012 shape functions. Figure 12 shows the 
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resulting shape. Although the designed shape did not fall 
right on top of the target shape, the residuals of the state 
and costate equations and the gradient of the cost function 
reached machine zero, which indicates that a minimum was 
reached. The cost function reached a value of 6 x 10 -9 . 




Fig. 11. Computed C p distribution for FX 60-126/1. 



Fig. 12. Test case 2. 

Next, an experiment was done to see how well the 
FX 60-126/1 airfoil can be represented with the NACA 0012 
shape functions. Figure 13 shows the result. The NACA 
0012 shape functions clearly do a good job everywhere except 
near the trailing edge. The reason that the optimum shape 
in the previous experiment does not correspond to the shape 
obtained from the shape fitting is not clear; one reason may be 
the poor quality of the grid because the airfoil is not smooth. 


FX60— 126/1 


Fitted with NACA 0012 shape functions 



Fig. 13. Shape fitting with 
NACA 0012 shape functions. 


Test Case 3. 

A third test was done; in this case, the fitted airfoil was 
used to generate the target potential. This shape is very close 
to the FX 60-126/1 airfoil and is smooth because it is based 
on smooth shape functions. The result of the design is shown 
in fig. 14. As expected, the final shape fell on top of the 
target shape. The residuals of the state and costate equations 
and the gradient of the cost function are shown in fig. 15. 



Target shape (Based on NACA 0012 shape functions) 

Final shape 



Fig. 14. Test case 3. 



Fig. 15. Convergence history. 
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The Efficiency of One-Shot Method 

Finally, the performance of the One-Shot method with 
respect to pure analysis is presented. The efficiency of a 
design method is defined as the ratio of the central processing 
unit (CPU) time that is required for the complete design 
process / n to the CPU time that is required to do a single 
analysis / i . Figure 16 shows this ratio to ft a plotted against 
the number of grid cells for the last test case. The figure 
shows that as the grid become finer the cost of design drops in 
comparison with the cost of one analysis. For the finest grid 
considered here, this ratio dropped below 4. The efficiencies 
for the other cases were similar. 



n 


Fig. 16. Efficiency of the One-Shot method. 
V. Concluding Remarks 


An efficient method for the design of optimal airfoil 
shapes has been presented in this paper. This method brings 
the cost of the design process to the same order as that of the 
analysis problem. It offers significant potential in the design 
of optimal aircraft configurations at a reasonable computer 
cost. However, much work is still required before practical 
use can be made of this method. 
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